#######################################################
#####REPLICATION ARCHIVE FOR NUANCED ACCOUNTABILITY####
#########BRITISH JOURNAL OF POLITICAL SCIENCE##########
########## INDIVIDUAL-LEVEL ANALYSES ##################
###de Kadt & Lieberman, 2017

#######################################################
################ SETUP BEGINS HERE ####################
#######################################################
## USER: Input your own working directory for files on line 50 ##
###Packages
library("foreign")
library("stargazer")
library("lmtest")
library("Zelig")
library("MatchIt")
library("rms")
library("ggplot2")
library("texreg")
library("car")
library("VGAM")
library("effects")
library("sandwich")
library("Matching")
library("ebal")
library("xtable")

###Cluster Robust SEs###
clse.f <- function(dat,fm, cluster){
  require(sandwich)
  require(lmtest)
  not <- attr(fm$model,"na.action")
  if( ! is.null(not)){
    cluster <- cluster[-not]
    dat <- dat[-not,]
  }
  with(dat,{
    M <- length(unique(cluster))
    N <- length(cluster)
    K <- fm$rank
    dfc <- (M/(M-1))*((N-1)/(N-K))
    uj <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
    vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
    coeftest(fm, vcovCL)
  }
  )
}

###Working Directories### 
dir = paste("####INPUT WORKING FILE DIRECTORY BETWEEN QUOTATION MARKS ####") 
setwd(paste(dir))

## Sets location of output
stgz <- paste(dir, "/output/", sep="")
fig <- paste(dir, "/output/", sep="")

afrob <- read.dta("nuanced_accountability_southafrica.dta") #keep original dataset
afrob$wave <- as.factor(afrob$wave)
afrob$catbcode <- as.factor(afrob$catbcode) #municipality code
afrob$w_log_pop_ward <- log(afrob$w_pop_ward)

##this variable is a 3 point scale: do you have toilet? water?
afrob$owntoil_water <- afrob$owntoilet + afrob$ownwater

##Removes missing categories
afrob$catbcode <- as.numeric(afrob$catbcode)
afrob$catbcode <- as.factor(afrob$catbcode)

afrob_ba <- subset(afrob, black == 1) #only blacks
afrob_b <- subset(afrob, black == 1 & seatshare >= .5) #only blacks in  ANC-controlled munis
afrob_bo <- subset(afrob, black == 1 & seatshare < .5) #only blacks in non ANC-controlled munis
afrob_b_ww <- subset(afrob_b, ownwater == 1) #subset of those blacks in ANC areas with water in house

#for analyses using aggregate data, use only years that are proximate to afrobarometer survey
afrob_b_agg <- subset(afrob_b, wave==906 | wave==903 | wave==902)
afrob_b_aggo <- subset(afrob_bo, wave==906 | wave==903 | wave==902)

## controls and controls for ward-level analyses
controls <- " + urban + female + cashjob + educ + age + home_zulu + home_xhosa" 
controls_w <- "districtcode + w_log_pop_ward + urban + female + cashjob + educ + age + home_zulu + home_xhosa" 


#### MAIN ANALYSES OF SOUTH AFRICA DATA BY SERVICE PROVISION SECTOR

### Coeff estimate OF effects of WATER
dvarnames <- c("govt_watsanit", "ownliving", "loc_coun_rate", "prez_rate", "anc_vote")
dvarlabels <- c("rate water/san", "rate own living", "rate local council", "rate president", "vote anc")

## create matrices to store the results
modelfits.ownwat <- vector(length(dvarnames), mode = "list") # matrix
names(modelfits.ownwat) <- dvarnames #name of each model based on dv name
modelfits.ownwat.c <- modelfits.ownwat # matrix to store clustered standard errors
names(modelfits.ownwat.c) <- dvarnames #again, name of each model

modelfits.eawat <- vector(length(dvarnames), mode = "list")
names(modelfits.eawat) <- dvarnames
modelfits.eawat.c <- modelfits.eawat
names(modelfits.eawat.c) <- dvarnames

modelfits.catwat <- vector(length(dvarnames), mode = "list")
names(modelfits.catwat) <- dvarnames
modelfits.catwat.c <- modelfits.catwat
names(modelfits.catwat.c) <- dvarnames

modelfits.wardwat <- vector(length(dvarnames), mode = "list")
names(modelfits.wardwat) <- dvarnames
modelfits.wardwat.c <- modelfits.wardwat
names(modelfits.wardwat.c) <- dvarnames

## Loop to go through each DV for "ownwater" (from AB), enumeration water (from AB), 
## catb water (cat_b or catbcode is for each municipality... ) and ward-level data.

for(i in dvarnames) { 
  
  modelformula.ownwat <- paste(i, " ~ ownwater + wave + catbcode + ",controls)
  modelfits.ownwat[[i]] <- lm(as.formula(modelformula.ownwat), data = afrob_b)
  modelfits.ownwat.c[[i]] <- clse.f(afrob_b, modelfits.ownwat[[i]], afrob_b$cat_b) 
  
  modelformula.eawat <- paste(i, " ~ ea_water + wave + catbcode + ",controls)
  modelfits.eawat[[i]] <- lm(as.formula(modelformula.eawat), data = afrob_b)
  modelfits.eawat.c[[i]] <- clse.f(afrob_b, modelfits.eawat[[i]], afrob_b$cat_b) 
  
  modelformula.catwat <- paste(i, " ~ home_water_nonwhite + wave  + ",controls)
  modelfits.catwat[[i]] <- lm(as.formula(modelformula.catwat), data = afrob_b_agg)
  modelfits.catwat.c[[i]] <- clse.f(afrob_b_agg, modelfits.catwat[[i]], afrob_b_agg$cat_b) 
  
  modelformula.wardwat <- paste(i, " ~ w_home_water_ward_nonwhite + wave  + ",controls_w)
  modelfits.wardwat[[i]] <- lm(as.formula(modelformula.wardwat), data = afrob_b_agg)
  modelfits.wardwat.c[[i]] <- clse.f(afrob_b_agg, modelfits.wardwat[[i]], afrob_b_agg$w_ward_id) 
}


### Coeff estimate OF REFUSE
##Basically same as above, but no afrobarometer data (own or ea) on refuse collection

dvarnames <- c("loc_govt_refuse", "ownliving", "loc_coun_rate", "prez_rate", "anc_vote")
dvarlabels <- c("rate \nrefuse", "rate own living", "rate local council", "rate president", "vote anc")


modelfits.catref <- vector(length(dvarnames), mode = "list")
names(modelfits.catref) <- dvarnames
modelfits.catref.c <- modelfits.catref
names(modelfits.catref.c) <- dvarnames

modelfits.wardref <- vector(length(dvarnames), mode = "list")
names(modelfits.wardref) <- dvarnames
modelfits.wardref.c <- modelfits.catref
names(modelfits.wardref.c) <- dvarnames


for(i in dvarnames) { 
  
  modelformula.catref <- paste(i, " ~ refuse_nonwhite  + ",controls)
  modelfits.catref[[i]] <- lm(as.formula(modelformula.catref), data = afrob_b_agg)
  modelfits.catref.c[[i]] <- clse.f(afrob_b_agg, modelfits.catref[[i]], afrob_b_agg$cat_b) 
  
  modelformula.wardref <- paste(i, " ~ w_refuse_ward_nonwhite   + ",controls_w)
  modelfits.wardref[[i]] <- lm(as.formula(modelformula.wardref), data = afrob_b_agg)
  modelfits.wardref.c[[i]] <- clse.f(afrob_b_agg, modelfits.wardref[[i]], afrob_b_agg$w_ward_id) 
}


### Coeff estimate OF SEWER/TOILET
dvarnames <- c("govt_watsanit", "ownliving", "loc_coun_rate", "prez_rate", "anc_vote")
dvarlabels <- c("rate water/san", "rate own living", "rate local council", "rate president", "vote anc")


modelfits.owntoi <- vector(length(dvarnames), mode = "list")
names(modelfits.owntoi) <- dvarnames
modelfits.owntoi.c <- modelfits.owntoi
names(modelfits.owntoi.c) <- dvarnames

modelfits.easew <- vector(length(dvarnames), mode = "list")
names(modelfits.easew) <- dvarnames
modelfits.easew.c <- modelfits.easew
names(modelfits.easew.c) <- dvarnames

modelfits.cattoi <- vector(length(dvarnames), mode = "list")
names(modelfits.cattoi) <- dvarnames
modelfits.cattoi.c <- modelfits.cattoi
names(modelfits.cattoi.c) <- dvarnames

modelfits.wardtoi <- vector(length(dvarnames), mode = "list")
names(modelfits.wardtoi) <- dvarnames
modelfits.wardtoi.c <- modelfits.wardtoi
names(modelfits.wardtoi.c) <- dvarnames


for(i in dvarnames) { 
  
  modelformula.owntoi <- paste(i, " ~ owntoilet + catbcode + ",controls)
  modelfits.owntoi[[i]] <- lm(as.formula(modelformula.owntoi), data = afrob_b)
  modelfits.owntoi.c[[i]] <- clse.f(afrob_b, modelfits.owntoi[[i]], afrob_b$cat_b) 
  
  modelformula.easew <- paste(i, " ~ ea_sewer + wave + catbcode + ",controls)
  modelfits.easew[[i]] <- lm(as.formula(modelformula.easew), data = afrob_b)
  modelfits.easew.c[[i]] <- clse.f(afrob_b, modelfits.easew[[i]], afrob_b$cat_b) 
  
  modelformula.cattoi <- paste(i, " ~ all_flush_nonwhite + wave  + ",controls)
  modelfits.cattoi[[i]] <- lm(as.formula(modelformula.cattoi), data = afrob_b_agg)
  modelfits.cattoi.c[[i]] <- clse.f(afrob_b_agg, modelfits.cattoi[[i]], afrob_b_agg$cat_b) 

  modelformula.wardtoi <- paste(i, " ~ w_all_flush_ward_nonwhite + wave  + ",controls_w)
  modelfits.wardtoi[[i]] <- lm(as.formula(modelformula.wardtoi), data = afrob_b_agg)
  modelfits.wardtoi.c[[i]] <- clse.f(afrob_b_agg, modelfits.wardtoi[[i]], afrob_b_agg$w_ward_id)   
}



#### REPEAT ANALYSES ABOVE... BUT FOR OPPOSITION-CONTROLLED MUNIS

### OPPOSITION AREAS ####
# -- repeat core analyses

#############################################################

### Coeff estimate OF effects of WATER
dvarnames <- c("govt_watsanit", "ownliving", "loc_coun_rate", "prez_rate", "anc_vote")
dvarlabels <- c("rate water/san", "rate own living", "rate local council", "rate president", "vote anc")

## create matrices to store the results
modelfits.ownwato <- vector(length(dvarnames), mode = "list") # matrix
names(modelfits.ownwato) <- dvarnames #name of each model based on dv name
modelfits.ownwato.c <- modelfits.ownwato # matrix to store clustered standard errors
names(modelfits.ownwato.c) <- dvarnames #again, name of each model

modelfits.eawato <- vector(length(dvarnames), mode = "list")
names(modelfits.eawato) <- dvarnames
modelfits.eawato.c <- modelfits.eawato
names(modelfits.eawato.c) <- dvarnames

modelfits.catwato <- vector(length(dvarnames), mode = "list")
names(modelfits.catwato) <- dvarnames
modelfits.catwato.c <- modelfits.catwato
names(modelfits.catwato.c) <- dvarnames

modelfits.wardwato <- vector(length(dvarnames), mode = "list")
names(modelfits.wardwato) <- dvarnames
modelfits.wardwato.c <- modelfits.wardwato
names(modelfits.wardwato.c) <- dvarnames

## Loop to go through each DV for "ownwater" (from AB), enumeration water (from AB), 
## catb water (cat_b or catbcode is for each municipality... ) and ward-level data.

for(i in dvarnames) { 
  
  modelformula.ownwato <- paste(i, " ~ ownwater + wave + catbcode + ",controls)
  modelfits.ownwato[[i]] <- lm(as.formula(modelformula.ownwato), data = afrob_bo)
  modelfits.ownwato.c[[i]] <- clse.f(afrob_bo, modelfits.ownwato[[i]], afrob_bo$cat_b) 
  
  modelformula.eawato <- paste(i, " ~ ea_water + wave + catbcode + ",controls)
  modelfits.eawato[[i]] <- lm(as.formula(modelformula.eawato), data = afrob_bo)
  modelfits.eawato.c[[i]] <- clse.f(afrob_bo, modelfits.eawato[[i]], afrob_bo$cat_b) 
  
  modelformula.catwato <- paste(i, " ~ home_water_nonwhite + wave  + ",controls)
  modelfits.catwato[[i]] <- lm(as.formula(modelformula.catwato), data = afrob_b_aggo)
  modelfits.catwato.c[[i]] <- clse.f(afrob_b_aggo, modelfits.catwato[[i]], afrob_b_aggo$cat_b) 
  
  modelformula.wardwato <- paste(i, " ~ w_home_water_ward_nonwhite + wave  + ",controls_w)
  modelfits.wardwato[[i]] <- lm(as.formula(modelformula.wardwato), data = afrob_b_aggo)
  modelfits.wardwato.c[[i]] <- clse.f(afrob_b_aggo, modelfits.wardwato[[i]], afrob_b_aggo$w_ward_id) 
}



### Coeff estimate OF REFUSE
##Basically same as above, but no afrobarometer data (own or ea) on refuse collection

dvarnames <- c("loc_govt_refuse", "ownliving", "loc_coun_rate", "prez_rate", "anc_vote")
dvarlabels <- c("rate \nrefuse", "rate own living", "rate local council", "rate president", "vote anc")


modelfits.catrefo <- vector(length(dvarnames), mode = "list")
names(modelfits.catrefo) <- dvarnames
modelfits.catrefo.c <- modelfits.catrefo
names(modelfits.catrefo.c) <- dvarnames

modelfits.wardrefo <- vector(length(dvarnames), mode = "list")
names(modelfits.wardrefo) <- dvarnames
modelfits.wardrefo.c <- modelfits.catrefo
names(modelfits.wardrefo.c) <- dvarnames


for(i in dvarnames) { 
  
  modelformula.catrefo <- paste(i, " ~ refuse_nonwhite  + ",controls)
  modelfits.catrefo[[i]] <- lm(as.formula(modelformula.catrefo), data = afrob_b_aggo)
  modelfits.catrefo.c[[i]] <- clse.f(afrob_b_aggo, modelfits.catrefo[[i]], afrob_b_aggo$cat_b) 
  
  modelformula.wardrefo <- paste(i, " ~ w_refuse_ward_nonwhite   + ",controls_w)
  modelfits.wardrefo[[i]] <- lm(as.formula(modelformula.wardrefo), data = afrob_b_aggo)
  modelfits.wardrefo.c[[i]] <- clse.f(afrob_b_aggo, modelfits.wardrefo[[i]], afrob_b_aggo$w_ward_id) 
}


### Coeff estimate OF SEWER/TOILET
dvarnames <- c("govt_watsanit", "ownliving", "loc_coun_rate", "prez_rate", "anc_vote")
dvarlabels <- c("rate water/san", "rate own living", "rate local council", "rate president", "vote anc")


modelfits.owntoio <- vector(length(dvarnames), mode = "list")
names(modelfits.owntoio) <- dvarnames
modelfits.owntoio.c <- modelfits.owntoio
names(modelfits.owntoio.c) <- dvarnames

modelfits.easewo <- vector(length(dvarnames), mode = "list")
names(modelfits.easewo) <- dvarnames
modelfits.easewo.c <- modelfits.easewo
names(modelfits.easewo.c) <- dvarnames

modelfits.cattoio <- vector(length(dvarnames), mode = "list")
names(modelfits.cattoio) <- dvarnames
modelfits.cattoio.c <- modelfits.cattoio
names(modelfits.cattoio.c) <- dvarnames

modelfits.wardtoio <- vector(length(dvarnames), mode = "list")
names(modelfits.wardtoio) <- dvarnames
modelfits.wardtoio.c <- modelfits.wardtoio
names(modelfits.wardtoio.c) <- dvarnames


for(i in dvarnames) { 
  
  modelformula.owntoio <- paste(i, " ~ owntoilet + catbcode + ",controls)
  modelfits.owntoio[[i]] <- lm(as.formula(modelformula.owntoio), data = afrob_bo)
  modelfits.owntoio.c[[i]] <- clse.f(afrob_bo, modelfits.owntoio[[i]], afrob_bo$cat_b) 
  
  modelformula.easewo <- paste(i, " ~ ea_sewer + wave + catbcode + ",controls)
  modelfits.easewo[[i]] <- lm(as.formula(modelformula.easewo), data = afrob_bo)
  modelfits.easewo.c[[i]] <- clse.f(afrob_bo, modelfits.easewo[[i]], afrob_bo$cat_b) 
  
  modelformula.cattoio <- paste(i, " ~ all_flush_nonwhite + wave  + ",controls)
  modelfits.cattoio[[i]] <- lm(as.formula(modelformula.cattoio), data = afrob_b_aggo)
  modelfits.cattoio.c[[i]] <- clse.f(afrob_b_aggo, modelfits.cattoio[[i]], afrob_b_aggo$cat_b) 
  
  modelformula.wardtoio <- paste(i, " ~ w_all_flush_ward_nonwhite + wave  + ",controls_w)
  modelfits.wardtoio[[i]] <- lm(as.formula(modelformula.wardtoio), data = afrob_b_aggo)
  modelfits.wardtoio.c[[i]] <- clse.f(afrob_b_aggo, modelfits.wardtoio[[i]], afrob_b_aggo$w_ward_id)   
}

####################################################
################ Paper Figure 2 ####################
####################################################
########## PLOTS OF ORIGINAL COEFFICIENT ESTIMATES
### ORGANIZED BY LEVEL OF ANALYSIS
pdf(file=paste(fig,"owneffab.pdf",sep=""), width = 4, height = 5)
par(mfrow = c(2,1), cex = .7)

dvarlabels <- c("rate \nwater/san", "rate \nown liv", "rate \nloc counc", "rate \nprez", "vote \nANC")
dvarnames <- c("govt_watsanit", "ownliving", "loc_coun_rate", "prez_rate", "anc_vote")


plot(c(1,2,3,4,5), c(modelfits.ownwat.c[[dvarnames[1]]][2,1],
                     modelfits.ownwat.c[[dvarnames[2]]][2,1],
                     modelfits.ownwat.c[[dvarnames[3]]][2,1],
                     modelfits.ownwat.c[[dvarnames[4]]][2,1],
                     modelfits.ownwat.c[[dvarnames[5]]][2,1]),
     ylim = c(-.3,.3),
     xlim= c(0,6),
     xaxt = "n",
     xlab= "",
     ylab = "",
     pch = 16,
     main = "Water (Household)\nANC-controlled munis"
)

axis(tick = FALSE, 1, 1:5, labels = dvarlabels, cex.axis=.8)

for (i in 1:5){
  lines(c(i,i), c(modelfits.ownwat.c[[dvarnames[i]]][2,1] - (1.96 * modelfits.ownwat.c[[dvarnames[i]]][2,2]),
                  modelfits.ownwat.c[[dvarnames[i]]][2,1] + (1.96 * modelfits.ownwat.c[[dvarnames[i]]][2,2])))
}
lines(c(0.0, 6.0), c(0,0), lty = 2)


plot(c(1,2,3,4,5), c(modelfits.owntoi.c[[dvarnames[1]]][2,1],
                     modelfits.owntoi.c[[dvarnames[2]]][2,1],
                     modelfits.owntoi.c[[dvarnames[3]]][2,1],
                     modelfits.owntoi.c[[dvarnames[4]]][2,1],
                     modelfits.owntoi.c[[dvarnames[5]]][2,1]),
     ylim = c(-.3,.3),
     
     xlim= c(0,6),
     xaxt = "n",
     xlab= "",
     ylab = "",
     pch = 16,
     main = "Toilet (Household)\nANC-controlled munis")

axis(tick = FALSE, 1, 1:5, labels = dvarlabels, cex.axis=.8)


for (i in 1:5){
  lines(c(i,i), c(modelfits.owntoi.c[[dvarnames[i]]][2,1] - (1.96 * modelfits.owntoi.c[[dvarnames[i]]][2,2]),
                  modelfits.owntoi.c[[dvarnames[i]]][2,1] + (1.96 * modelfits.owntoi.c[[dvarnames[i]]][2,2])))
}
lines(c(0.0, 6.0), c(0,0), lty = 2)

dev.off()


### EA LEVEL
pdf(file=paste(fig,"EAeffab.pdf",sep=""), width = 4, height = 5)
par(mfrow = c(2,1), cex = .7)


dvarnames <- c("govt_watsanit", "ownliving", "loc_coun_rate", "prez_rate", "anc_vote")
dvarlabels <- c("rate \nwater/san", "rate \nown liv", "rate \nloc counc", "rate \nprez", "vote \nANC")

plot(c(1,2,3,4,5), c(modelfits.eawat.c[[dvarnames[1]]][2,1],
                     modelfits.eawat.c[[dvarnames[2]]][2,1],
                     modelfits.eawat.c[[dvarnames[3]]][2,1],
                     modelfits.eawat.c[[dvarnames[4]]][2,1],
                     modelfits.eawat.c[[dvarnames[5]]][2,1]),
     ylim = c(-.3,.3),
     xlim= c(0,6),
     xaxt = "n",
     xlab= "",
     ylab = "",
     pch = 16,
     main = "Water (Enumeration area)\nANC-controlled munis")

axis(tick = FALSE, 1, 1:5, labels = dvarlabels, cex.axis=.8)


for (i in 1:5){
  lines(c(i,i), c(modelfits.eawat.c[[dvarnames[i]]][2,1] - (1.96 * modelfits.eawat.c[[dvarnames[i]]][2,2]),
                  modelfits.eawat.c[[dvarnames[i]]][2,1] + (1.96 * modelfits.eawat.c[[dvarnames[i]]][2,2])))
}
lines(c(0.0, 6.0), c(0,0), lty = 2)


plot(c(1,2,3,4,5), c(modelfits.easew.c[[dvarnames[1]]][2,1],
                     modelfits.easew.c[[dvarnames[2]]][2,1],
                     modelfits.easew.c[[dvarnames[3]]][2,1],
                     modelfits.easew.c[[dvarnames[4]]][2,1],
                     modelfits.easew.c[[dvarnames[5]]][2,1]),
     ylim = c(-.3,.3),
     
     xlim= c(0,6),
     xaxt = "n",
     xlab= "",
     ylab = "",
     pch = 16,
     main = "Sewer (Enumeration area)\nANC-controlled munis")

axis(tick = FALSE, 1, 1:5,  labels = dvarlabels, cex.axis=.8)


for (i in 1:5){
  lines(c(i,i), c(modelfits.easew.c[[dvarnames[i]]][2,1] - (1.96 * modelfits.easew.c[[dvarnames[i]]][2,2]),
                  modelfits.easew.c[[dvarnames[i]]][2,1] + (1.96 * modelfits.easew.c[[dvarnames[i]]][2,2])))
}
lines(c(0.0, 6.0), c(0,0), lty = 2)


dev.off()



########## FIGURE 2 CONTINUED -- PLOTS FOR OPPOSITION AREAS

pdf(file=paste(fig,"owneffabo.pdf",sep=""), width = 4, height = 5)
par(mfrow = c(2,1), cex = .7)

dvarlabels <- c("rate \nwater/san", "rate \nown liv", "rate \nloc counc", "rate \nprez", "vote \nANC")
dvarnames <- c("govt_watsanit", "ownliving", "loc_coun_rate", "prez_rate", "anc_vote")


plot(c(1,2,3,4,5), c(modelfits.ownwato.c[[dvarnames[1]]][2,1],
                     modelfits.ownwato.c[[dvarnames[2]]][2,1],
                     modelfits.ownwato.c[[dvarnames[3]]][2,1],
                     modelfits.ownwato.c[[dvarnames[4]]][2,1],
                     modelfits.ownwato.c[[dvarnames[5]]][2,1]),
     ylim = c(-.8,.3),
     xlim= c(0,6),
     xaxt = "n",
     xlab= "",
     ylab = "",
     pch = 16,
     main = "Water (Household)\nOpposition-controlled munis"
)

axis(tick = FALSE, 1, 1:5, labels = dvarlabels, cex.axis=.8)

for (i in 1:5){
  lines(c(i,i), c(modelfits.ownwato.c[[dvarnames[i]]][2,1] - (1.96 * modelfits.ownwato.c[[dvarnames[i]]][2,2]),
                  modelfits.ownwato.c[[dvarnames[i]]][2,1] + (1.96 * modelfits.ownwato.c[[dvarnames[i]]][2,2])))
}
lines(c(0.0, 6.0), c(0,0), lty = 2)


plot(c(1,2,3,4,5), c(modelfits.owntoio.c[[dvarnames[1]]][2,1],
                     modelfits.owntoio.c[[dvarnames[2]]][2,1],
                     modelfits.owntoio.c[[dvarnames[3]]][2,1],
                     modelfits.owntoio.c[[dvarnames[4]]][2,1],
                     modelfits.owntoio.c[[dvarnames[5]]][2,1]),
     ylim = c(-.8,.3),
     
     xlim= c(0,6),
     xaxt = "n",
     xlab= "",
     ylab = "",
     pch = 16,
     main = "Toilet (Household)\nOpposition-controlled munis")

axis(tick = FALSE, 1, 1:5, labels = dvarlabels, cex.axis=.8)


for (i in 1:5){
  lines(c(i,i), c(modelfits.owntoio.c[[dvarnames[i]]][2,1] - (1.96 * modelfits.owntoio.c[[dvarnames[i]]][2,2]),
                  modelfits.owntoio.c[[dvarnames[i]]][2,1] + (1.96 * modelfits.owntoio.c[[dvarnames[i]]][2,2])))
}
lines(c(0.0, 6.0), c(0,0), lty = 2)


dev.off()


### EA LEVEL
pdf(file=paste(fig,"EAeffabo.pdf",sep=""), width = 4, height = 5)
par(mfrow = c(2,1), cex = .7)


dvarnames <- c("govt_watsanit", "ownliving", "loc_coun_rate", "prez_rate", "anc_vote")
dvarlabels <- c("rate \nwater/san", "rate \nown liv", "rate \nloc counc", "rate \nprez", "vote \nANC")

plot(c(1,2,3,4,5), c(modelfits.eawato.c[[dvarnames[1]]][2,1],
                     modelfits.eawato.c[[dvarnames[2]]][2,1],
                     modelfits.eawato.c[[dvarnames[3]]][2,1],
                     modelfits.eawato.c[[dvarnames[4]]][2,1],
                     modelfits.eawato.c[[dvarnames[5]]][2,1]),
     ylim = c(-.6,.3),
     xlim= c(0,6),
     xaxt = "n",
     xlab= "",
     ylab = "",
     pch = 16,
     main = "Water (Enumeration area)\nOpposition-controlled munis")

axis(tick = FALSE, 1, 1:5, labels = dvarlabels, cex.axis=.8)


for (i in 1:5){
  lines(c(i,i), c(modelfits.eawato.c[[dvarnames[i]]][2,1] - (1.96 * modelfits.eawato.c[[dvarnames[i]]][2,2]),
                  modelfits.eawato.c[[dvarnames[i]]][2,1] + (1.96 * modelfits.eawato.c[[dvarnames[i]]][2,2])))
}
lines(c(0.0, 6.0), c(0,0), lty = 2)


plot(c(1,2,3,4,5), c(modelfits.easewo.c[[dvarnames[1]]][2,1],
                     modelfits.easewo.c[[dvarnames[2]]][2,1],
                     modelfits.easewo.c[[dvarnames[3]]][2,1],
                     modelfits.easewo.c[[dvarnames[4]]][2,1],
                     modelfits.easewo.c[[dvarnames[5]]][2,1]),
     ylim = c(-.3,.3),
     
     xlim= c(0,6),
     xaxt = "n",
     xlab= "",
     ylab = "",
     pch = 16,
     main = "Sewer (Enumeration area)\nOpposition-controlled munis")

axis(tick = FALSE, 1, 1:5,  labels = dvarlabels, cex.axis=.8)


for (i in 1:5){
  lines(c(i,i), c(modelfits.easewo.c[[dvarnames[i]]][2,1] - (1.96 * modelfits.easewo.c[[dvarnames[i]]][2,2]),
                  modelfits.easewo.c[[dvarnames[i]]][2,1] + (1.96 * modelfits.easewo.c[[dvarnames[i]]][2,2])))
}
lines(c(0.0, 6.0), c(0,0), lty = 2)


dev.off()

####################################################
###### NOTE USER: Matching / Figure 3 at end #######
####################################################

####################################################
################ Paper Figure 4 ####################
####################################################
### MUNI LEVEL
pdf(file=paste(fig,"MUNIeffab.pdf",sep=""), width = 5, height = 7)
par(mfrow = c(3,1), cex = .7)
dvarnames <- c("govt_watsanit", "ownliving", "loc_coun_rate", "prez_rate", "anc_vote")
dvarlabels <- c("rate \nwater/san", "rate \nown liv", "rate \nloc counc", "rate \nprez", "vote \nANC")

plot(c(1,2,3,4,5), c(modelfits.catwat.c[[dvarnames[1]]][2,1],
                     modelfits.catwat.c[[dvarnames[2]]][2,1],
                     modelfits.catwat.c[[dvarnames[3]]][2,1],
                     modelfits.catwat.c[[dvarnames[4]]][2,1],
                     modelfits.catwat.c[[dvarnames[5]]][2,1]),
     ylim = c(-.5,.5),
     xlim= c(0,6),
     xaxt = "n",
     xlab= "",
     ylab = "",
     pch = 16,
     main = "Water  (Municipal coverage)")

axis(tick = FALSE, 1, 1:5, labels = dvarlabels, cex.axis=.8)


for (i in 1:5){
  lines(c(i,i), c(modelfits.catwat.c[[dvarnames[i]]][2,1] - (1.96 * modelfits.catwat.c[[dvarnames[i]]][2,2]),
                  modelfits.catwat.c[[dvarnames[i]]][2,1] + (1.96 * modelfits.catwat.c[[dvarnames[i]]][2,2])))
}
lines(c(0.0, 6.0), c(0,0), lty = 2)



plot(c(1,2,3,4,5), c(modelfits.cattoi.c[[dvarnames[1]]][2,1],
                     modelfits.cattoi.c[[dvarnames[2]]][2,1],
                     modelfits.cattoi.c[[dvarnames[3]]][2,1],
                     modelfits.cattoi.c[[dvarnames[4]]][2,1],
                     modelfits.cattoi.c[[dvarnames[5]]][2,1]),
     ylim = c(-.3,.3),
     xlim= c(0,6),
     xaxt = "n",
     xlab= "",
     ylab = "",
     pch = 16,
     main = "Toilet (Municipal coverage)")

axis(tick = FALSE, 1, 1:5, labels = dvarlabels, cex.axis=.8)

for (i in 1:5){
  lines(c(i,i), c(modelfits.cattoi.c[[dvarnames[i]]][2,1] - (1.96 * modelfits.cattoi.c[[dvarnames[i]]][2,2]),
                  modelfits.cattoi.c[[dvarnames[i]]][2,1] + (1.96 * modelfits.cattoi.c[[dvarnames[i]]][2,2])))
}
lines(c(0.0, 6.0), c(0,0), lty = 2)


dvarnames <- c("loc_govt_refuse", "ownliving", "loc_coun_rate", "prez_rate", "anc_vote")
dvarlabels <- c("rate \nrefuse", "rate \nown liv", "rate \nloc counc", "rate \nprez", "vote \nANC")



plot(c(1,2,3,4,5), c(modelfits.catref.c[[dvarnames[1]]][2,1],
                     modelfits.catref.c[[dvarnames[2]]][2,1],
                     modelfits.catref.c[[dvarnames[3]]][2,1],
                     modelfits.catref.c[[dvarnames[4]]][2,1],
                     modelfits.catref.c[[dvarnames[5]]][2,1]),
     ylim = c(-.3,.3),
     xlim= c(0,6),
     xaxt = "n",
     xlab= "",
     ylab = "",
     pch = 16,
     main = "Refuse (Municipal coverage)")

axis(tick = FALSE, 1, 1:5, labels = dvarlabels, cex.axis=.8)

for (i in 1:5){
  lines(c(i,i), c(modelfits.catref.c[[dvarnames[i]]][2,1] - (1.96 * modelfits.catref.c[[dvarnames[i]]][2,2]),
                  modelfits.catref.c[[dvarnames[i]]][2,1] + (1.96 * modelfits.catref.c[[dvarnames[i]]][2,2])))
}
lines(c(0.0, 6.0), c(0,0), lty = 2)

dev.off()

### WARD LEVEL
pdf(file=paste(fig,"WARDeffab.pdf",sep=""), width = 5, height = 7)
par(mfrow = c(3,1), cex = .7)
dvarnames <- c("govt_watsanit", "ownliving", "loc_coun_rate", "prez_rate", "anc_vote")
dvarlabels <- c("rate \nwater/san", "rate \nown liv", "rate \nloc counc", "rate \nprez", "vote \nANC")

plot(c(1,2,3,4,5), c(modelfits.wardwat.c[[dvarnames[1]]][2,1],
                     modelfits.wardwat.c[[dvarnames[2]]][2,1],
                     modelfits.wardwat.c[[dvarnames[3]]][2,1],
                     modelfits.wardwat.c[[dvarnames[4]]][2,1],
                     modelfits.wardwat.c[[dvarnames[5]]][2,1]),
     ylim = c(-.3,.3),
     xlim= c(0,6),
     xaxt = "n",
     xlab= "",
     ylab = "",
     pch = 16,
     main = "Water  (Ward coverage)")

axis(tick = FALSE, 1, 1:5, labels = dvarlabels, cex.axis=.8)


for (i in 1:5){
  lines(c(i,i), c(modelfits.wardwat.c[[dvarnames[i]]][2,1] - (1.96 * modelfits.wardwat.c[[dvarnames[i]]][2,2]),
                  modelfits.wardwat.c[[dvarnames[i]]][2,1] + (1.96 * modelfits.wardwat.c[[dvarnames[i]]][2,2])))
}
lines(c(0.0, 6.0), c(0,0), lty = 2)



plot(c(1,2,3,4,5), c(modelfits.wardtoi.c[[dvarnames[1]]][2,1],
                     modelfits.wardtoi.c[[dvarnames[2]]][2,1],
                     modelfits.wardtoi.c[[dvarnames[3]]][2,1],
                     modelfits.wardtoi.c[[dvarnames[4]]][2,1],
                     modelfits.wardtoi.c[[dvarnames[5]]][2,1]),
     ylim = c(-.3,.3),
     xlim= c(0,6),
     xaxt = "n",
     xlab= "",
     ylab = "",
     pch = 16,
     main = "Toilet (Ward coverage)")

axis(tick = FALSE, 1, 1:5, labels = dvarlabels, cex.axis=.8)

for (i in 1:5){
  lines(c(i,i), c(modelfits.wardtoi.c[[dvarnames[i]]][2,1] - (1.96 * modelfits.wardtoi.c[[dvarnames[i]]][2,2]),
                  modelfits.wardtoi.c[[dvarnames[i]]][2,1] + (1.96 * modelfits.wardtoi.c[[dvarnames[i]]][2,2])))
}
lines(c(0.0, 6.0), c(0,0), lty = 2)


dvarnames <- c("loc_govt_refuse", "ownliving", "loc_coun_rate", "prez_rate", "anc_vote")
dvarlabels <- c("rate \nrefuse", "rate \nown liv", "rate \nloc counc", "rate \nprez", "vote \nANC")


plot(c(1,2,3,4,5), c(modelfits.wardref.c[[dvarnames[1]]][2,1],
                     modelfits.wardref.c[[dvarnames[2]]][2,1],
                     modelfits.wardref.c[[dvarnames[3]]][2,1],
                     modelfits.wardref.c[[dvarnames[4]]][2,1],
                     modelfits.wardref.c[[dvarnames[5]]][2,1]),
     ylim = c(-.3,.3),
     xlim= c(0,6),
     xaxt = "n",
     xlab= "",
     ylab = "",
     pch = 16,
     main = "Refuse (Ward coverage)")

axis(tick = FALSE, 1, 1:5, labels = dvarlabels, cex.axis=.8)

for (i in 1:5){
  lines(c(i,i), c(modelfits.wardref.c[[dvarnames[i]]][2,1] - (1.96 * modelfits.wardref.c[[dvarnames[i]]][2,2]),
                  modelfits.wardref.c[[dvarnames[i]]][2,1] + (1.96 * modelfits.wardref.c[[dvarnames[i]]][2,2])))
}
lines(c(0.0, 6.0), c(0,0), lty = 2)


dev.off()

####################################################
################# CORRUPTION ANALYSIS ##############
###### Paper Table 6 (plus appendix sec. 8.3) ######
####################################################
afrob_b.906 = subset(afrob_b, wave==906) #only final wave has both toilet and water data

corrupt1 <- lm(anc_vote ~ corrupt + urban + female + cashjob + educ + age,
               data = afrob_b.906)
corrupt1.c <- clse.f(afrob_b.906, corrupt1, afrob_b.906$cat_b) 

corrupt2 <- lm(corrupt ~ ownwater + urban + female + cashjob + educ + age,
               data = afrob_b.906)
corrupt2.c <- clse.f(afrob_b.906, corrupt2, afrob_b.906$cat_b) 

corrupt3 <- lm(corrupt ~ owntoilet + urban + female + cashjob + educ + age,
               data = afrob_b.906)
corrupt3.c <- clse.f(afrob_b.906, corrupt3, afrob_b.906$cat_b) 

corrupt.list <- list(corrupt1.c, corrupt2.c, corrupt3.c)

corrupt.list.clse = lapply(corrupt.list, `[`,,2)


tblab <- paste(stgz,"corrupt.tex", sep="") #ADJUST LABEL
stargazer(corrupt1, corrupt2, corrupt3,
          style = "ajps",
          #         font.size = "tiny",
          dep.var.labels.include = FALSE,
          column.labels = c("ANC Vote", "Perceive Corruption", "Perceive Corruption"),
          omit.stat = c("f", "adj.rsq", "ser"),
          title = "Service Provision and Perceptions of Corruption (AB Round 5 data, black respondents in ANC-controlled municipalities)",
          notes ="Clustered (municipality) standard errors in parentheses.",
          se = (c(corrupt.list.clse[1], corrupt.list.clse[2], corrupt.list.clse[3])),
          out = tblab,
          label= "corrupt"
)


####################################################
########## RELATIVE DEPRIVATION ANALYSIS ###########
########### (Appendix Section 8.4) #################
####################################################
afrob_b_agg_nw = subset(afrob_b_agg, ownwater==0)
reldepw <- lm(anc_vote ~ home_water + whitefrac + urban + female + cashjob + educ +
                age, data = afrob_b_agg_nw)
reldepw.c <- clse.f(afrob_b_agg_nw, reldepw, afrob_b_agg_nw$cat_b) 

reldepw2 <- lm(ownliving ~ home_water + whitefrac + urban + female + cashjob + educ +
                 age , data = afrob_b_agg_nw)
reldepw2.c <- clse.f(afrob_b_agg_nw, reldepw2, afrob_b_agg_nw$cat_b) 

reldepw3 <- lm(liv_compoth ~ home_water + whitefrac + urban + female + cashjob + educ +
                 age , data = afrob_b_agg_nw)
reldepw3.c <- clse.f(afrob_b_agg_nw, reldepw3, afrob_b_agg_nw$cat_b) 


afrob_b_agg_nt = subset(afrob_b_agg, owntoilet<.50)
reldept <- lm(anc_vote ~ all_flush + whitefrac + urban + female + cashjob + educ +
                age , data = afrob_b_agg_nt)
reldept.c <- clse.f(afrob_b_agg_nt, reldept, afrob_b_agg_nt$cat_b) 

reldept2 <- lm(ownliving ~ all_flush + whitefrac + urban + female + cashjob + educ +
                 age , data = afrob_b_agg_nt)
reldept2.c <- clse.f(afrob_b_agg_nt, reldept2, afrob_b_agg_nt$cat_b) 

reldept3 <- lm(liv_compoth ~ all_flush + whitefrac + urban + female + cashjob + educ +
                 age , data = afrob_b_agg_nt)
reldept3.c <- clse.f(afrob_b_agg_nt, reldept3, afrob_b_agg_nt$cat_b) 


# for reporting clustered SE's
reldep.list <- list(reldepw2.c, reldept2.c, reldepw3.c, reldept3.c, reldepw.c, reldept.c)

reldep.list.clse = lapply(reldep.list, `[`,,2)

tblab <- paste(stgz,"reldep.tex", sep="") #ADJUST LABEL
stargazer(reldepw2, reldept2, reldepw3, reldept3, reldepw, reldept,
          style = "ajps",
          font.size = "tiny",
          dep.var.labels.include = FALSE,
          column.labels = c("Own Living", "Own Living", "Liv Comp Othrs", "Liv Comp Othrs","ANC Vote",  "ANC Vote"),
          omit.stat = c("f", "adj.rsq", "ser"),
          title = "The effects of municipal-level service delivery coverage among those without household services",
          notes ="Clustered (municipality) standard errors in parentheses.",
          se = c(reldep.list.clse[1], reldep.list.clse[2], reldep.list.clse[3], reldep.list.clse[4], 
                 reldep.list.clse[5], reldep.list.clse[6]),
          out = tblab,
          label = "reldep"
)


####################################################
############## DISAPPOINTMENT ANALYSIS #############
#### Is outcome driven by quality of services? #####
########### (Appendix Section 8.5) #################
####################################################
### Coeff estimate OF WATER QUALITY
dvarnames <- c("govt_watsanit", "ownliving", "loc_coun_rate", "prez_rate", "anc_vote")
dvarlabels <- c("rate water/san", "rate own living", "rate local council", "rate president", "vote anc")

modelfits.qwat <- vector(length(dvarnames), mode = "list")
names(modelfits.qwat) <- dvarnames
modelfits.qwat.c <- modelfits.qwat
names(modelfits.qwat.c) <- dvarnames

for(i in dvarnames) { 
  modelformula.qwat <- paste(i, " ~ enoughwater + wave + catbcode + ",controls)
  modelfits.qwat[[i]] <- lm(as.formula(modelformula.qwat), data = afrob_b_ww)
  modelfits.qwat.c[[i]] <- clse.f(afrob_b, modelfits.qwat[[i]], afrob_b_ww$cat_b) 
  
}

##WATER QUALITY
pdf(file=paste(stgz, "waterqualab.pdf", sep=""), width = 5, height = 5)
par(mfrow = c(1,1), cex = .6)

dvarnames <- c("govt_watsanit", "ownliving", "loc_coun_rate", "prez_rate", "anc_vote")
dvarlabels <- c("rate \nwater/san", "rate \nown liv", "rate \nloc counc", "rate \nprez", "vote \nANC")

plot(c(1,2,3,4,5), c(modelfits.qwat.c[[dvarnames[1]]][2,1],
                     modelfits.qwat.c[[dvarnames[2]]][2,1],
                     modelfits.qwat.c[[dvarnames[3]]][2,1],
                     modelfits.qwat.c[[dvarnames[4]]][2,1],
                     modelfits.qwat.c[[dvarnames[5]]][2,1]),
     ylim = c(-.3,.3), pch = 16,
     xlim= c(0,6),
     xaxt = "n",
     xlab= "",
     ylab = "Coeff estimate of quality of water supply (among those with water)"
)

axis(tick = FALSE, 1, 1:5, labels = dvarlabels, cex.axis=.75)

for (i in 1:5){
  lines(c(i,i), c(modelfits.qwat.c[[dvarnames[i]]][2,1] - (1.96 * modelfits.qwat.c[[dvarnames[i]]][2,2]),
                  modelfits.qwat.c[[dvarnames[i]]][2,1] + (1.96 * modelfits.qwat.c[[dvarnames[i]]][2,2])))
}
lines(c(0.0, 6.0), c(0,0), lty = 2)

dev.off()

## DOES HAVING WATER AFFECT LIKELIHOOD OF PRIORITIZING WATER? (Reported in main text)
modwatpri <- lm(waterpriority ~ ownwater + urban + female + cashjob + educ +
                   age + as.factor(wave), data = afrob_b)

modwatpri.c <- clse.f(afrob_b, modwatpri, afrob_b$cat_b) 

wptab <- table(afrob_b$ownwater, afrob_b$waterpriority)
wptab

#Water as a priority amongst those without water:
wptab[1,2] / (wptab[1,1] + wptab[1,2])

#Water as a priority amongst those with water:
(wptab[2,2] + wptab[3,2] ) / (wptab[2,1] + wptab[2,2] + wptab[3,1] + wptab[3,2])

## Who prioritizes water, electric, basic services?
basictab <- table(afrob_ba$servpri)
basictab[2] / (basictab[1] + basictab[2])

# And if in low service area?
basictabl <- table(afrob_ba$servpri[afrob_ba$ea_water==0 & afrob_ba$ea_elec==0])
basictabl[2] / (basictabl[1] + basictabl[2])

#################################
###### Tables For Appendix ######
#################################
## For reporting clustered SE's
ownwat.cl.se.out = lapply(modelfits.ownwat.c, `[`,,2)
eawat.cl.se.out = lapply(modelfits.eawat.c, `[`,,2)
catwat.cl.se.out = lapply(modelfits.catwat.c, `[`,,2)
wardwat.cl.se.out = lapply(modelfits.wardwat.c, `[`,,2)
catref.cl.se.out = lapply(modelfits.catref.c, `[`,,2)
wardref.cl.se.out = lapply(modelfits.wardref.c, `[`,,2)
owntoi.cl.se.out = lapply(modelfits.owntoi.c, `[`,,2)
easew.cl.se.out = lapply(modelfits.easew.c, `[`,,2)
cattoi.cl.se.out = lapply(modelfits.cattoi.c, `[`,,2)
wardtoi.cl.se.out = lapply(modelfits.wardtoi.c, `[`,,2)

#create table label for stargazer output
tablen <- 0

#Models
mnames <- c("ownwat", "eawat", "catwat", "wardwat", "catref", "wardref", "owntoi", "easew", "cattoi", "wardtoi")
tblabels <- c("Household water access", "Enumeration Area water access", "Municipal Water Coverage (Percent Households)", "Ward Water Coverage (Percent Households)",
              "Municipal Refuse Collection Coverage (Percent Households)", "Ward Refuse Collection Coverage (Percent Households)",
              "Household toilet access", "Enumeration Area sewer line access", "Municipal Toilet Service Coverage (Percent Households)", "Ward Toilet Service Coverage (Percent Households)")

tablen <- tablen + 1
tblab <- paste(stgz,"table", tablen, "append.tex", sep="") #ADJUST LABEL

stargazer(modelfits.ownwat,
          style = "ajps",
          font.size = "tiny",
          dep.var.labels.include = FALSE,
          column.labels = c("Rate service", "Own living", "Rate local gov", "Rate prez", "Vote ANC"),
          omit.stat = c("f", "adj.rsq", "ser"),
          title = paste("OLS estimates of effect of ", tblabels[tablen]), 
          se = ownwat.cl.se.out,
          omit = "catbcode",
          notes = "Clustered (muni) standard errors in parentheses. ANC-controlled areas",
          out = tblab
)

tablen <- tablen + 1
tblab <- paste(stgz,"table", tablen, "append.tex", sep="") #ADJUST LABEL

stargazer(modelfits.eawat,
          style = "ajps",
          font.size = "tiny",
          dep.var.labels.include = FALSE,
          column.labels = c("Rate service", "Own living", "Rate local gov", "Rate prez", "Vote ANC"),
          omit.stat = c("f", "adj.rsq", "ser"),
          title = paste("OLS estimates of effect of ", tblabels[tablen]), 
          se = eawat.cl.se.out,
          omit = "catbcode",
          notes = "Clustered (muni) standard errors in parentheses. ANC-controlled areas",
          out = tblab
)

tablen <- tablen + 1
tblab <- paste(stgz,"table", tablen, "append.tex", sep="") #ADJUST LABEL

stargazer(modelfits.catwat,
          style = "ajps",
          font.size = "tiny",
          dep.var.labels.include = FALSE,
          column.labels = c("Rate service", "Own living", "Rate local gov", "Rate prez", "Vote ANC"),
          omit.stat = c("f", "adj.rsq", "ser"),
          title = paste("OLS estimates of effect of ", tblabels[tablen]), 
          se = catwat.cl.se.out,
          omit = "catbcode",
          notes = "Clustered (muni) standard errors in parentheses. ANC-controlled areas",
          out = tblab
)
tablen <- tablen + 1
tblab <- paste(stgz,"table", tablen, "append.tex", sep="") #ADJUST LABEL

stargazer(modelfits.wardwat,
          style = "ajps",
          font.size = "tiny",
          dep.var.labels.include = FALSE,
          column.labels = c("Rate service", "Own living", "Rate local gov", "Rate prez", "Vote ANC"),
          omit.stat = c("f", "adj.rsq", "ser"),
          title = paste("OLS estimates of effect of ", tblabels[tablen]), 
          se = wardwat.cl.se.out,
          omit = "catbcode",
          notes = "Clustered (muni) standard errors in parentheses. ANC-controlled areas",
          out = tblab
)
tablen <- tablen + 1
tblab <- paste(stgz,"table", tablen, "append.tex", sep="") #ADJUST LABEL

stargazer(modelfits.catref,
          style = "ajps",
          font.size = "tiny",
          dep.var.labels.include = FALSE,
          column.labels = c("Rate service", "Own living", "Rate local gov", "Rate prez", "Vote ANC"),
          omit.stat = c("f", "adj.rsq", "ser"),
          title = paste("OLS estimates of effect of ", tblabels[tablen]), 
          se = catref.cl.se.out,
          omit = "catbcode",
          notes = "Clustered (muni) standard errors in parentheses. ANC-controlled areas",
          out = tblab
)
tablen <- tablen + 1
tblab <- paste(stgz,"table", tablen, "append.tex", sep="") #ADJUST LABEL

stargazer(modelfits.wardref,
          style = "ajps",
          font.size = "tiny",
          dep.var.labels.include = FALSE,
          column.labels = c("Rate service", "Own living", "Rate local gov", "Rate prez", "Vote ANC"),
          omit.stat = c("f", "adj.rsq", "ser"),
          title = paste("OLS estimates of effect of ", tblabels[tablen]), 
          se = wardref.cl.se.out,
          omit = "catbcode",
          notes = "Clustered (muni) standard errors in parentheses. ANC-controlled areas",
          out = tblab
)

tablen <- tablen + 1
tblab <- paste(stgz,"table", tablen, "append.tex", sep="") #ADJUST LABEL

stargazer(modelfits.owntoi,
          style = "ajps",
          font.size = "tiny",
          dep.var.labels.include = FALSE,
          column.labels = c("Rate service", "Own living", "Rate local gov", "Rate prez", "Vote ANC"),
          omit.stat = c("f", "adj.rsq", "ser"),
          title = paste("OLS estimates of effect of ", tblabels[tablen]), 
          se = owntoi.cl.se.out,
          omit = "catbcode",
          notes = "Clustered (muni) standard errors in parentheses. ANC-controlled areas",
          out = tblab
)
tablen <- tablen + 1
tblab <- paste(stgz,"table", tablen, "append.tex", sep="") #ADJUST LABEL

stargazer(modelfits.easew,
          style = "ajps",
          font.size = "tiny",
          dep.var.labels.include = FALSE,
          column.labels = c("Rate service", "Own living", "Rate local gov", "Rate prez", "Vote ANC"),
          omit.stat = c("f", "adj.rsq", "ser"),
          title = paste("OLS estimates of effect of ", tblabels[tablen]), 
          se = easew.cl.se.out,
          omit = "catbcode",
          notes = "Clustered (muni) standard errors in parentheses. ANC-controlled areas",
          out = tblab
)
tablen <- tablen + 1
tblab <- paste(stgz,"table", tablen, "append.tex", sep="") #ADJUST LABEL

stargazer(modelfits.cattoi,
          style = "ajps",
          font.size = "tiny",
          dep.var.labels.include = FALSE,
          column.labels = c("Own Living", "Own Living", "Liv Comp Othrs", "Liv Comp Othrs","ANC Vote",  "ANC Vote"),
          omit.stat = c("f", "adj.rsq", "ser"),
          title = paste("OLS estimates of effect of ", tblabels[tablen]), 
          se = cattoi.cl.se.out,
          omit = "catbcode",
          notes = "Clustered (muni) standard errors in parentheses. ANC-controlled areas",
          out = tblab
)
tablen <- tablen + 1
tblab <- paste(stgz,"table", tablen, "append.tex", sep="") #ADJUST LABEL

stargazer(modelfits.wardtoi,
          style = "ajps",
          font.size = "tiny",
          dep.var.labels.include = FALSE,
          column.labels = c("Rate service", "Own living", "Rate local gov", "Rate prez", "Vote ANC"),
          omit.stat = c("f", "adj.rsq", "ser"),
          title = paste("OLS estimates of effect of ", tblabels[tablen]), 
          se = wardtoi.cl.se.out,
          omit = "catbcode",
          notes = "Clustered (muni) standard errors in parentheses. ANC-controlled areas",
          out = tblab
)

#################################
## Tables For Opposition Areas ##
#################################
## For reporting clustered SE's
ownwato.cl.se.out = lapply(modelfits.ownwato.c, `[`,,2)
eawato.cl.se.out = lapply(modelfits.eawato.c, `[`,,2)
catwato.cl.se.out = lapply(modelfits.catwato.c, `[`,,2)
wardwato.cl.se.out = lapply(modelfits.wardwato.c, `[`,,2)
catrefo.cl.se.out = lapply(modelfits.catrefo.c, `[`,,2)
wardrefo.cl.se.out = lapply(modelfits.wardrefo.c, `[`,,2)
owntoio.cl.se.out = lapply(modelfits.owntoio.c, `[`,,2)
easewo.cl.se.out = lapply(modelfits.easewo.c, `[`,,2)
cattoio.cl.se.out = lapply(modelfits.cattoio.c, `[`,,2)
wardtoio.cl.se.out = lapply(modelfits.wardtoio.c, `[`,,2)

#create table label for stargazer output
tablen <- 0

#Opposition-area models
mnames <- c("ownwato", "eawato", "catwato", "wardwato", "catrefo", "wardrefo", "owntoio", "easewo", "cattoio", "wardtoio")
tblabels <- c("Household water access", "Enumeration Area water access", "Municipal Water Coverage (Percent Households)", "Ward Water Coverage (Percent Households)",
              "Municipal Refuse Collection Coverage (Percent Households)", "Ward Refuse Collection Coverage (Percent Households)",
              "Household toilet access", "Enumeration Area sewer line access", "Municipal Toilet Service Coverage (Percent Households)", "Ward Toilet Service Coverage (Percent Households)")

tablen <- tablen + 1
tblab <- paste(stgz,"tableo", tablen, "append.tex", sep="") #ADJUST LABEL

stargazer(modelfits.ownwato,
          style = "ajps",
          font.size = "tiny",
          dep.var.labels.include = FALSE,
          column.labels = c("Rate service", "Own living", "Rate local gov", "Rate prez", "Vote ANC"),
          omit.stat = c("f", "adj.rsq", "ser"),
          title = paste("OLS estimates of effect of ", tblabels[tablen]), 
          se = ownwato.cl.se.out,
          omit = "catbcode",
          notes = "Clustered (muni) standard errors in parentheses. Non-ANC areas",
          out = tblab
)

tablen <- tablen + 1
tblab <- paste(stgz,"tableo", tablen, "append.tex", sep="") #ADJUST LABEL

stargazer(modelfits.eawato,
          style = "ajps",
          font.size = "tiny",
          dep.var.labels.include = FALSE,
          column.labels = c("Rate service", "Own living", "Rate local gov", "Rate prez", "Vote ANC"),
          omit.stat = c("f", "adj.rsq", "ser"),
          title = paste("OLS estimates of effect of ", tblabels[tablen]), 
          se = eawato.cl.se.out,
          omit = "catbcode",
          notes = "Clustered (muni) standard errors in parentheses. Non-ANC areas",
          out = tblab
)

tablen <- tablen + 1
tblab <- paste(stgz,"tableo", tablen, "append.tex", sep="") #ADJUST LABEL

stargazer(modelfits.catwato,
          style = "ajps",
          font.size = "tiny",
          dep.var.labels.include = FALSE,
          column.labels = c("Rate service", "Own living", "Rate local gov", "Rate prez", "Vote ANC"),
          omit.stat = c("f", "adj.rsq", "ser"),
          title = paste("OLS estimates of effect of ", tblabels[tablen]), 
          se = catwato.cl.se.out,
          omit = "catbcode",
          notes = "Clustered (muni) standard errors in parentheses. Non-ANC areas",
          out = tblab
)
tablen <- tablen + 1
tblab <- paste(stgz,"tableo", tablen, "append.tex", sep="") #ADJUST LABEL

stargazer(modelfits.wardwato,
          style = "ajps",
          font.size = "tiny",
          dep.var.labels.include = FALSE,
          column.labels = c("Rate service", "Own living", "Rate local gov", "Rate prez", "Vote ANC"),
          omit.stat = c("f", "adj.rsq", "ser"),
          title = paste("OLS estimates of effect of ", tblabels[tablen]), 
          se = wardwato.cl.se.out,
          omit = "catbcode",
          notes = "Clustered (muni) standard errors in parentheses. Non-ANC areas",
          out = tblab
)
tablen <- tablen + 1
tblab <- paste(stgz,"tableo", tablen, "append.tex", sep="") #ADJUST LABEL

stargazer(modelfits.catrefo,
          style = "ajps",
          font.size = "tiny",
          dep.var.labels.include = FALSE,
          column.labels = c("Rate service", "Own living", "Rate local gov", "Rate prez", "Vote ANC"),
          omit.stat = c("f", "adj.rsq", "ser"),
          title = paste("OLS estimates of effect of ", tblabels[tablen]), 
          se = catrefo.cl.se.out,
          omit = "catbcode",
          notes = "Clustered (muni) standard errors in parentheses. Non-ANC areas",
          out = tblab
)
tablen <- tablen + 1
tblab <- paste(stgz,"tableo", tablen, "append.tex", sep="") #ADJUST LABEL

stargazer(modelfits.wardrefo,
          style = "ajps",
          font.size = "tiny",
          dep.var.labels.include = FALSE,
          column.labels = c("Rate service", "Own living", "Rate local gov", "Rate prez", "Vote ANC"),
          omit.stat = c("f", "adj.rsq", "ser"),
          title = paste("OLS estimates of effect of ", tblabels[tablen]), 
          se = wardrefo.cl.se.out,
          omit = "catbcode",
          notes = "Clustered (muni) standard errors in parentheses. Non-ANC areas",
          out = tblab
)

tablen <- tablen + 1
tblab <- paste(stgz,"tableo", tablen, "append.tex", sep="") #ADJUST LABEL

stargazer(modelfits.owntoio,
          style = "ajps",
          font.size = "tiny",
          dep.var.labels.include = FALSE,
          column.labels = c("Rate service", "Own living", "Rate local gov", "Rate prez", "Vote ANC"),
          omit.stat = c("f", "adj.rsq", "ser"),
          title = paste("OLS estimates of effect of ", tblabels[tablen]), 
          se = owntoio.cl.se.out,
          omit = "catbcode",
          notes = "Clustered (muni) standard errors in parentheses. Non-ANC areas",
          out = tblab
)
tablen <- tablen + 1
tblab <- paste(stgz,"tableo", tablen, "append.tex", sep="") #ADJUST LABEL

stargazer(modelfits.easewo,
          style = "ajps",
          font.size = "tiny",
          dep.var.labels.include = FALSE,
          column.labels = c("Rate service", "Own living", "Rate local gov", "Rate prez", "Vote ANC"),
          omit.stat = c("f", "adj.rsq", "ser"),
          title = paste("OLS estimates of effect of ", tblabels[tablen]), 
          se = easewo.cl.se.out,
          omit = "catbcode",
          notes = "Clustered (muni) standard errors in parentheses. Non-ANC areas",
          out = tblab
)
tablen <- tablen + 1
tblab <- paste(stgz,"tableo", tablen, "append.tex", sep="") #ADJUST LABEL

stargazer(modelfits.cattoio,
          style = "ajps",
          font.size = "tiny",
          dep.var.labels.include = FALSE,
          column.labels = c("Own Living", "Own Living", "Liv Comp Othrs", "Liv Comp Othrs","ANC Vote",  "ANC Vote"),
          omit.stat = c("f", "adj.rsq", "ser"),
          title = paste("OLS estimates of effect of ", tblabels[tablen]), 
          se = cattoio.cl.se.out,
          omit = "catbcode",
          notes = "Clustered (muni) standard errors in parentheses. Non-ANC areas",
          out = tblab
)
tablen <- tablen + 1
tblab <- paste(stgz,"tableo", tablen, "append.tex", sep="") #ADJUST LABEL

stargazer(modelfits.wardtoio,
          style = "ajps",
          font.size = "tiny",
          dep.var.labels.include = FALSE,
          column.labels = c("Rate service", "Own living", "Rate local gov", "Rate prez", "Vote ANC"),
          omit.stat = c("f", "adj.rsq", "ser"),
          title = paste("OLS estimates of effect of ", tblabels[tablen]), 
          se = wardtoio.cl.se.out,
          omit = "catbcode",
          notes = "Clustered (muni) standard errors in parentheses. Non-ANC areas",
          out = tblab
)

###############################################
####### Botswana, Lesotho, and Namibia ########
##### Paper Table 5 and aAppendix Sec 10 ######
###############################################
######BOTSWANA
afrob_b <- read.dta("nuanced_accountability_botswana.dta") #keep original dataset
afrob_b$round <- as.factor(afrob_b$round)

#create table number for stargazer output
tablen <- 0

## controls 
controls <- " + urban + female + cashjob + educ + age + home_setswana" 

#####ANALYSES#####
### Coeff estimate OF effects of WATER
dvarnames <- c("govt_watsanit", "ownliving", "loc_coun_rate", "prez_rate", "bdp_vote")
dvarlabels <- c("rate water/san", "rate own living", "rate local council", "rate president", "vote bdp")

## create matrices to store the results
modelfits.ownwat <- vector(length(dvarnames), mode = "list") # matrix
names(modelfits.ownwat) <- dvarnames #name of each model based on dv name

modelfits.eawat <- vector(length(dvarnames), mode = "list")
names(modelfits.eawat) <- dvarnames

## Loop to go through each DV for "ownwater" (from AB), enumeration water (from AB), 
for(i in dvarnames) { 
  modelformula.ownwat <- paste(i, " ~ ownwater + round  + ",controls)
  modelfits.ownwat[[i]] <- lm(as.formula(modelformula.ownwat), data = afrob_b)

  modelformula.eawat <- paste(i, " ~ ea_water + round  + ",controls)
  modelfits.eawat[[i]] <- lm(as.formula(modelformula.eawat), data = afrob_b)
}

tablen <- tablen + 1
tblab <- paste(stgz,"table", tablen, "bots.tex", sep="") #ADJUST LABEL
stargazer(modelfits.ownwat,
          style = "ajps",
          font.size = "tiny",
          dep.var.labels.include = FALSE,
          column.labels = dvarlabels,
          omit.stat = c("f", "adj.rsq", "ser"),
          title = "Botswana: Estimates of Household-level access to water",
          notes ="Afrobarometer data",
          out = tblab
)

tablen <- tablen + 1
tblab <- paste(stgz,"table", tablen, "bots.tex", sep="") #ADJUST LABEL
stargazer(modelfits.eawat,
          style = "ajps",
          font.size = "tiny",
          dep.var.labels.include = FALSE,
          column.labels = dvarlabels,
          omit.stat = c("f", "adj.rsq", "ser"),
          title = "Botswana: Estimates of Enumeration area access to water",
          notes ="Afrobarometer data",
          out = tblab
)

### Coeff estimate OF SEWER/TOILET
dvarnames <- c("govt_watsanit", "ownliving", "loc_coun_rate", "prez_rate", "bdp_vote")
dvarlabels <- c("rate water/san", "rate own living", "rate local council", "rate president", "vote bdp")

modelfits.owntoi <- vector(length(dvarnames), mode = "list")
names(modelfits.owntoi) <- dvarnames

modelfitseasew <- vector(length(dvarnames), mode = "list")
names(modelfitseasew) <- dvarnames

for(i in dvarnames) { 
  modelformula.owntoi <- paste(i, " ~ owntoilet + ",controls)
  modelfits.owntoi[[i]] <- lm(as.formula(modelformula.owntoi), data = afrob_b)
  
  modelformulaeasew <- paste(i, " ~ ea_sewer + round + ",controls)
  modelfitseasew[[i]] <- lm(as.formula(modelformulaeasew), data = afrob_b)  
}

tablen <- tablen + 1
tblab <- paste(stgz,"table", tablen, "bots.tex", sep="") #ADJUST LABEL
stargazer(modelfits.owntoi,
          style = "ajps",
          font.size = "tiny",
          dep.var.labels.include = FALSE,
          column.labels = dvarlabels,
          omit.stat = c("f", "adj.rsq", "ser"),
          title = "Botswana: Estimates of Household-level access to flush toilet",
          notes ="Afrobarometer data",
          out = tblab
)

tablen <- tablen + 1
tblab <- paste(stgz,"table", tablen, "bots.tex", sep="") #ADJUST LABEL
stargazer(modelfitseasew,
          style = "ajps",
          font.size = "tiny",
          dep.var.labels.include = FALSE,
          column.labels = dvarlabels,
          omit.stat = c("f", "adj.rsq", "ser"),
          title = "Botswana: Estimates of Enumeration area sewer line",
          notes ="Afrobarometer data",
          out = tblab
)

######LESOTHO
afrob_b <- read.dta("nuanced_accountability_lesotho.dta") #get original dataset
afrob_b$round <- as.factor(afrob_b$round)

## create table number for stargazer output
tablen <- 0

## controls 
controls <- " + urban + female + cashjob + educ + age" 

#####ANALYSES#####
### Coeff estimate OF effects of WATER
dvarnames <- c("govt_watsanit", "ownliving", "loc_coun_rate", "prez_rate", "lcd_vote")
dvarlabels <- c("rate water/san", "rate own living", "rate local council", "rate president", "vote lcd")

## create matrices to store the results
modelfits.eawat <- vector(length(dvarnames), mode = "list")
names(modelfits.eawat) <- dvarnames

## Loop to go through each DV for "ownwater" (from AB), enumeration water (from AB), 
for(i in dvarnames) {   
  modelformula.eawat <- paste(i, " ~ ea_water  + ",controls)
  modelfits.eawat[[i]] <- lm(as.formula(modelformula.eawat), data = afrob_b)
}

tablen <- tablen + 1
tblab <- paste(stgz,"table", tablen, "les.tex", sep="") #ADJUST LABEL
stargazer(modelfits.eawat,
          style = "ajps",
          font.size = "tiny",
          dep.var.labels.include = FALSE,
          column.labels = dvarlabels,
          omit.stat = c("f", "adj.rsq", "ser"),
          title = "Lesotho: Estimates of Enumeration area access to water",
          notes ="Afrobarometer data",
          out = tblab
)

### Coeff estimate OF SEWER/TOILET
dvarnames <- c("govt_watsanit", "ownliving", "loc_coun_rate", "prez_rate", "lcd_vote")
dvarlabels <- c("rate water/san", "rate own living", "rate local council", "rate president", "vote lcd")

modelfitseasew <- vector(length(dvarnames), mode = "list")
names(modelfitseasew) <- dvarnames

for(i in dvarnames) { 
  modelformulaeasew <- paste(i, " ~ ea_sewer + ",controls)
  modelfitseasew[[i]] <- lm(as.formula(modelformulaeasew), data = afrob_b)
}

tablen <- tablen + 1
tblab <- paste(stgz,"table", tablen, "les.tex", sep="") #ADJUST LABEL
stargazer(modelfitseasew,
          style = "ajps",
          font.size = "tiny",
          dep.var.labels.include = FALSE,
          column.labels = dvarlabels,
          omit.stat = c("f", "adj.rsq", "ser"),
          title = "Lesotho: Estimates of Enumeration area sewer line",
          notes ="Afrobarometer data",
          out = tblab
)


######NAMIBIA
afrob_b <- read.dta("nuanced_accountability_namibia.dta") #get original dataset
afrob_b$round <- as.factor(afrob_b$round)

## create table number for stargazer output
tablen <- 0

## controls 
controls <- " + urban + female + cashjob + educ + age + home_ovambo" 

#####ANALYSES#####
### Coeff estimate OF effects of WATER
dvarnames <- c("govt_watsanit", "ownliving", "loc_coun_rate", "prez_rate", "swapo_vote")
dvarlabels <- c("rate water/san", "rate own living", "rate local council", "rate president", "vote swapo")

## create matrices to store the results
modelfits.ownwat <- vector(length(dvarnames), mode = "list") # matrix
names(modelfits.ownwat) <- dvarnames #name of each model based on dv name

modelfits.eawat <- vector(length(dvarnames), mode = "list")
names(modelfits.eawat) <- dvarnames

## Loop to go through each DV for "ownwater" (from AB), enumeration water (from AB), 
for(i in dvarnames) {   
  modelformula.ownwat <- paste(i, " ~ ownwater + round  + ",controls)
  modelfits.ownwat[[i]] <-lm(as.formula(modelformula.ownwat), data = afrob_b)
  
  modelformula.eawat <- paste(i, " ~ ea_water + round  + ",controls)
  modelfits.eawat[[i]] <- lm(as.formula(modelformula.eawat), data = afrob_b)
}

tablen <- tablen + 1
tblab <- paste(stgz,"table", tablen, "nam.tex", sep="") #ADJUST LABEL
stargazer(modelfits.ownwat,
          style = "ajps",
          font.size = "tiny",
          dep.var.labels.include = FALSE,
          column.labels = dvarlabels,
          omit.stat = c("f", "adj.rsq", "ser"),
          title = "Namibia: Estimates of Household-level access to water",
          notes ="Afrobarometer data",
          out = tblab
)

tablen <- tablen + 1
tblab <- paste(stgz,"table", tablen, "nam.tex", sep="") #ADJUST LABEL
stargazer(modelfits.eawat,
          style = "ajps",
          font.size = "tiny",
          dep.var.labels.include = FALSE,
          column.labels = dvarlabels,
          omit.stat = c("f", "adj.rsq", "ser"),
          title = "Namibia: Estimates of Enumeration area access to water",
          notes ="Afrobarometer data",
          out = tblab
)

### Coeff estimate OF SEWER/TOILET
dvarnames <- c("govt_watsanit", "ownliving", "loc_coun_rate", "prez_rate", "swapo_vote")
dvarlabels <- c("rate water/san", "rate own living", "rate local council", "rate president", "vote swapo")

modelfits.owntoi <- vector(length(dvarnames), mode = "list")
names(modelfits.owntoi) <- dvarnames

modelfitseasew <- vector(length(dvarnames), mode = "list")
names(modelfitseasew) <- dvarnames

for(i in dvarnames) {   
  modelformula.owntoi <- paste(i, " ~ owntoilet + ",controls)
  modelfits.owntoi[[i]] <- lm(as.formula(modelformula.owntoi),  data = afrob_b)

  modelformulaeasew <- paste(i, " ~ ea_sewer + round + ",controls)
  modelfitseasew[[i]] <- lm(as.formula(modelformulaeasew),  data = afrob_b)
}

tablen <- tablen + 1
tblab <- paste(stgz,"table", tablen, "nam.tex", sep="") #ADJUST LABEL
stargazer(modelfits.owntoi,
          style = "ajps",
          font.size = "tiny",
          dep.var.labels.include = FALSE,
          column.labels = dvarlabels,
          omit.stat = c("f", "adj.rsq", "ser"),
          title = "Namibia: Estimates of Household-level access to flush toilet",
          notes ="Afrobarometer data",
          out = tblab
)

tablen <- tablen + 1
tblab <- paste(stgz,"table", tablen, "nam.tex", sep="") #ADJUST LABEL
stargazer(modelfitseasew,
          style = "ajps",
          font.size = "tiny",
          dep.var.labels.include = FALSE,
          column.labels = dvarlabels,
          omit.stat = c("f", "adj.rsq", "ser"),
          title = "Namibia: Estimates of Enumeration area sewer line",
          notes ="Afrobarometer data",
          out = tblab
)

####################################################
############### Matching Analyses ##################
###### Paper Figure 3 and Appendix Section 9 #######
####################################################
##USER: Note, need to generate dataset separately for each sector so we don't lose observations from prior matching

###### WATER
modelfits.w <- matrix(data = NA, nrow = 5, ncol = 4, byrow = FALSE,
                      dimnames = NULL)

afrob <- read.dta("nuanced_accountability_southafrica.dta") #keep original dataset
afrob_b <- subset(afrob, black == 1 & seatshare > .5 & (wave == 905 | wave == 906))

afrob_b$cat_b_num = as.numeric(factor(afrob_b$cat_b))
afrob_b$district = as.numeric(as.factor(afrob_b$district))

store <- c('anc_vote', 'cashjob', 'home_eng', 'home_afrik', 'home_xhosa', 'home_zulu', 'ownwater_t', 'wave', 'age', 'log_pop', 'educ', 'educ_r', 'district', 'province', 'cat_b_num', 
           'urban', 'female', 'govt_watsanit', 'ownliving', 'loc_coun_rate', 'prez_rate')
afrob_match <- afrob_b[, store]
afrob_match_na <- na.omit(afrob_match)
afrob_b <- afrob_match_na

dvarnames <- c("govt_watsanit", "ownliving", "loc_coun_rate", "prez_rate", "anc_vote") #####Cs
dvarlabels <- c("rate water/san", "rate own living", "rate local council", "rate president", "vote anc")

# Estimate the ATE as the difference in means of ownwater
mean(afrob_b$govt_watsanit[afrob_b$ownwater_t == 1], na.rm=TRUE) - mean(afrob_b$govt_watsanit[afrob_b$ownwater_t == 0], na.rm=TRUE)
# And estimate its standard error
sqrt(var(afrob_b$govt_watsanit[afrob_b$ownwater_t == 1], na.rm=TRUE)/length(afrob_b$govt_watsanit[afrob_b$ownwater_t == 1]) + 
       var(afrob_b$govt_watsanit[afrob_b$ownwater_t == 0], na.rm=TRUE)/length(afrob_b$govt_watsanit[afrob_b$ownwater_t == 0]))
# Check balance in unmatched data
bal.formula = formula(ownwater_t ~ log_pop + age + educ_r + cashjob + urban + female + home_afrik + home_xhosa + home_zulu, na.rm = TRUE)
# Print a balance table
mb.unmatched = MatchBalance(bal.formula, data = afrob_b, print.level=0)
tab.unmatched = baltest.collect(mb.unmatched, var.names = c('log_pop', 'age', 'educ_r', 'cashjob', 'urban', 'female', "home_afrik", "home_xhosa", "home_zulu"), after=F)
print(xtable(tab.unmatched[, 1:6],  label='tab:unmatched-bal',caption='Covariate Balance in Unmatched Data'), caption.placement='top')


# Match on a set of covariates
vars = c('age', 'cat_b_num', 'log_pop', 'educ_r', 'cashjob', 'home_afrik', 'home_xhosa', 'home_zulu')
# Adjust for bias; function defaults will estimate the ATT with M=1
m.first = Match(Y=afrob_b$govt_watsanit, Tr=afrob_b$ownwater_t, X = afrob_b[, vars], BiasAdjust=TRUE)
summary(m.first)# treatment effect
# Check balance
mb.matched = MatchBalance(bal.formula, data=afrob_b, match.out=m.first, print.level=0)
# ATT and its SE can be found in the est and se elements of m.first
# Print a balance table
tab.matched = baltest.collect(mb.matched, var.names = c('log_pop', 'age', 'educ_r', 'cashjob', 'urban', 'female', "home_afrik", "home_xhosa", "home_zulu"))
print(xtable(tab.matched[, 1:7], caption='Covariate Balance in Matched Data', label='tab:matched-bal'), caption.placement='top')


# Match exactly on the geo IDs and cashjob and education
vars_exact = c('cat_b_num', 'urban', 'cashjob', 'educ_r', 'wave', "home_zulu" )

m.exact = Match(Y=afrob_b$govt_watsanit, Tr=afrob_b$ownwater_t, X = afrob_b[, vars_exact], BiasAdjust=TRUE, exact=T)
summary(m.exact) ## Treatment effect
# Check balance
mb.matched.exact = MatchBalance(bal.formula, data=afrob_b, match.out=m.exact, print.level=0)
# ATT and its SE can be found in the est and se elements of m.exact
# Print the table
tab.exact = baltest.collect(mb.matched.exact, var.names = c('log_pop', 'age', 'educ_r', 'cashjob', 'urban', 'female', "home_afrik", "home_xhosa", "home_zulu") )
print(xtable(tab.exact[, 1:7], caption='Covariate Balance in Exactly Matched Data', label='tab:exact-bal'), caption.placement='top')
modelfits.w[1,] <- c(m.exact$est, m.exact$se, m.exact$orig.nobs, m.exact$wnobs)

m.exact = Match(Y=afrob_b$ownliving, Tr=afrob_b$ownwater_t, X = afrob_b[, vars_exact], BiasAdjust=TRUE, exact=T)
summary(m.exact) ## Treatment effect
# Check balance
mb.matched.exact = MatchBalance(bal.formula, data=afrob_b, match.out=m.exact, print.level=0)
# ATT and its SE can be found in the est and se elements of m.exact
# Print the table
tab.exact = baltest.collect(mb.matched.exact, var.names = c('log_pop', 'age', 'educ_r', 'cashjob', 'urban', 'female', "home_afrik", "home_xhosa", "home_zulu") )
print(xtable(tab.exact[, 1:7], caption='Covariate Balance in Exactly Matched Data', label='tab:exact-bal'), caption.placement='top')
modelfits.w[2,] <- c(m.exact$est, m.exact$se, m.exact$orig.nobs, m.exact$wnobs)


m.exact = Match(Y=afrob_b$loc_coun_rate, Tr=afrob_b$ownwater_t, X = afrob_b[, vars_exact], BiasAdjust=TRUE, exact=T)
summary(m.exact) ## Treatment effect
# Check balance
mb.matched.exact = MatchBalance(bal.formula, data=afrob_b, match.out=m.exact, print.level=0)
# ATT and its SE can be found in the est and se elements of m.exact
# Print the table
tab.exact = baltest.collect(mb.matched.exact, var.names = c('log_pop', 'age', 'educ_r', 'cashjob', 'urban', 'female', "home_afrik", "home_xhosa", "home_zulu") )
print(xtable(tab.exact[, 1:7], caption='Covariate Balance in Exactly Matched Data', label='tab:exact-bal'), caption.placement='top')
modelfits.w[3,] <- c(m.exact$est, m.exact$se, m.exact$orig.nobs, m.exact$wnobs)


m.exact = Match(Y=afrob_b$prez_rate, Tr=afrob_b$ownwater_t, X = afrob_b[, vars_exact], BiasAdjust=TRUE, exact=T)
summary(m.exact) ## Treatment effect
# Check balance
mb.matched.exact = MatchBalance(bal.formula, data=afrob_b, match.out=m.exact, print.level=0)
# ATT and its SE can be found in the est and se elements of m.exact
# Print the table
tab.exact = baltest.collect(mb.matched.exact, var.names = c('log_pop', 'age', 'educ_r', 'cashjob', 'urban', 'female', "home_afrik", "home_xhosa", "home_zulu") )
print(xtable(tab.exact[, 1:7], caption='Covariate Balance in Exactly Matched Data', label='tab:exact-bal'), caption.placement='top')
modelfits.w[4,] <- c(m.exact$est, m.exact$se, m.exact$orig.nobs, m.exact$wnobs)

m.exact = Match(Y=afrob_b$anc_vote, Tr=afrob_b$ownwater_t, X = afrob_b[, vars_exact], BiasAdjust=TRUE, exact=T)
summary(m.exact) ## Treatment effect
# Check balance
mb.matched.exact = MatchBalance(bal.formula, data=afrob_b, match.out=m.exact, print.level=0)
# ATT and its SE can be found in the est and se elements of m.exact
# Print the table
tab.exact = baltest.collect(mb.matched.exact, var.names = c('log_pop', 'age', 'educ_r', 'cashjob', 'urban', 'female', "home_afrik", "home_xhosa", "home_zulu") )
print(xtable(tab.exact[, 1:7], caption='Covariate Balance in Exactly Matched Data', label='tab:exact-bal'), caption.placement='top')
modelfits.w[5,] <- c(m.exact$est, m.exact$se, m.exact$orig.nobs, m.exact$wnobs)


###### TOILET
modelfits.t <- matrix(data = NA, nrow = 5, ncol = 4, byrow = FALSE,
                      dimnames = NULL)

afrob <- read.dta("nuanced_accountability_southafrica.dta") #keep original dataset
afrob_b <- subset(afrob, black == 1 & seatshare > .5 & wave == 906)

afrob_b$cat_b_num = as.numeric(factor(afrob_b$cat_b))
afrob_b$district = as.numeric(as.factor(afrob_b$district))

store <- c('anc_vote', 'cashjob', 'home_eng', 'home_afrik', 'home_xhosa', 'home_zulu', 'owntoilet_b', 'wave', 'age', 'log_pop', 'educ', 'educ_r', 'district', 'province', 'cat_b_num', 
           'urban', 'female', 'govt_watsanit', 'ownliving', 'loc_coun_rate', 'prez_rate')
afrob_match <- afrob_b[, store]
afrob_match_na <- na.omit(afrob_match)
afrob_b <- afrob_match_na

dvarnames <- c("govt_watsanit", "ownliving", "loc_coun_rate", "prez_rate", "anc_vote") ###Cs
dvarlabels <- c("rate water/san", "rate own living", "rate local council", "rate president", "vote anc")


# Estimate the ATE as the difference in means of ownwater
mean(afrob_b$govt_watsanit[afrob_b$owntoilet_b == 1], na.rm=TRUE) - mean(afrob_b$govt_watsanit[afrob_b$owntoilet_b == 0], na.rm=TRUE)
# And estimate its standard error
sqrt(var(afrob_b$govt_watsanit[afrob_b$owntoilet_b == 1], na.rm=TRUE)/length(afrob_b$govt_watsanit[afrob_b$owntoilet_b == 1]) + 
       var(afrob_b$govt_watsanit[afrob_b$owntoilet_b == 0], na.rm=TRUE)/length(afrob_b$govt_watsanit[afrob_b$owntoilet_b == 0]))
# Check balance in unmatched data
bal.formula = formula(owntoilet_b ~ log_pop + age + educ_r + cashjob + urban + female + home_afrik + home_xhosa + home_zulu, na.rm = TRUE)
# Print a balance table
mb.unmatched = MatchBalance(bal.formula, data = afrob_b, print.level=0)
tab.unmatched = baltest.collect(mb.unmatched, var.names = c('log_pop', 'age', 'educ_r', 'cashjob', 'urban', 'female', "home_afrik", "home_xhosa", "home_zulu"), after=F)
print(xtable(tab.unmatched[, 1:6],  label='tab:unmatched-bal',caption='Covariate Balance in Unmatched Data'), caption.placement='top')


# Match on a set of covariates
vars = c('age', 'cat_b_num', 'log_pop', 'educ_r', 'cashjob', 'home_afrik', 'home_xhosa', 'home_zulu')
# Adjust for bias; function defaults will estimate the ATT with M=1
m.first = Match(Y=afrob_b$govt_watsanit, Tr=afrob_b$owntoilet_b, X = afrob_b[, vars], BiasAdjust=TRUE)
summary(m.first)# treatment effect
# Check balance
mb.matched = MatchBalance(bal.formula, data=afrob_b, match.out=m.first, print.level=0)
# ATT and its SE can be found in the est and se elements of m.first
# Print a balance table
tab.matched = baltest.collect(mb.matched, var.names = c('log_pop', 'age', 'educ_r', 'cashjob', 'urban', 'female', "home_afrik", "home_xhosa", "home_zulu"))
print(xtable(tab.matched[, 1:7], caption='Covariate Balance in Matched Data', label='tab:matched-bal'), caption.placement='top')


m.exact = Match(Y=afrob_b$govt_watsanit, Tr=afrob_b$owntoilet_b, X = afrob_b[, vars_exact], BiasAdjust=TRUE, exact=T)
summary(m.exact) ## Treatment effect
# Check balance
mb.matched.exact = MatchBalance(bal.formula, data=afrob_b, match.out=m.exact, print.level=0)
# ATT and its SE can be found in the est and se elements of m.exact
# Print the table
tab.exact = baltest.collect(mb.matched.exact, var.names = c('log_pop', 'age', 'educ_r', 'cashjob', 'urban', 'female', "home_afrik", "home_xhosa", "home_zulu") )
print(xtable(tab.exact[, 1:7], caption='Covariate Balance in Exactly Matched Data', label='tab:exact-bal'), caption.placement='top')
modelfits.t[1,] <- c(m.exact$est, m.exact$se, m.exact$orig.nobs, m.exact$wnobs)

m.exact = Match(Y=afrob_b$ownliving, Tr=afrob_b$owntoilet_b, X = afrob_b[, vars_exact], BiasAdjust=TRUE, exact=T)
summary(m.exact) ## Treatment effect
# Check balance
mb.matched.exact = MatchBalance(bal.formula, data=afrob_b, match.out=m.exact, print.level=0)
# ATT and its SE can be found in the est and se elements of m.exact
# Print the table
tab.exact = baltest.collect(mb.matched.exact, var.names = c('log_pop', 'age', 'educ_r', 'cashjob', 'urban', 'female', "home_afrik", "home_xhosa", "home_zulu") )
print(xtable(tab.exact[, 1:7], caption='Covariate Balance in Exactly Matched Data', label='tab:exact-bal'), caption.placement='top')
modelfits.t[2,] <- c(m.exact$est, m.exact$se, m.exact$orig.nobs, m.exact$wnobs)


m.exact = Match(Y=afrob_b$loc_coun_rate, Tr=afrob_b$owntoilet_b, X = afrob_b[, vars_exact], BiasAdjust=TRUE, exact=T)
summary(m.exact) ## Treatment effect
# Check balance
mb.matched.exact = MatchBalance(bal.formula, data=afrob_b, match.out=m.exact, print.level=0)
# ATT and its SE can be found in the est and se elements of m.exact
# Print the table
tab.exact = baltest.collect(mb.matched.exact, var.names = c('log_pop', 'age', 'educ_r', 'cashjob', 'urban', 'female', "home_afrik", "home_xhosa", "home_zulu") )
print(xtable(tab.exact[, 1:7], caption='Covariate Balance in Exactly Matched Data', label='tab:exact-bal'), caption.placement='top')
modelfits.t[3,] <- c(m.exact$est, m.exact$se, m.exact$orig.nobs, m.exact$wnobs)


m.exact = Match(Y=afrob_b$prez_rate, Tr=afrob_b$owntoilet_b, X = afrob_b[, vars_exact], BiasAdjust=TRUE, exact=T)
summary(m.exact) ## Treatment effect
# Check balance
mb.matched.exact = MatchBalance(bal.formula, data=afrob_b, match.out=m.exact, print.level=0)
# ATT and its SE can be found in the est and se elements of m.exact
# Print the table
tab.exact = baltest.collect(mb.matched.exact, var.names = c('log_pop', 'age', 'educ_r', 'cashjob', 'urban', 'female', "home_afrik", "home_xhosa", "home_zulu") )
print(xtable(tab.exact[, 1:7], caption='Covariate Balance in Exactly Matched Data', label='tab:exact-bal'), caption.placement='top')
modelfits.t[4,] <- c(m.exact$est, m.exact$se, m.exact$orig.nobs, m.exact$wnobs)

m.exact = Match(Y=afrob_b$anc_vote, Tr=afrob_b$owntoilet_b, X = afrob_b[, vars_exact], BiasAdjust=TRUE, exact=T)
summary(m.exact) ## Treatment effect
# Check balance
mb.matched.exact = MatchBalance(bal.formula, data=afrob_b, match.out=m.exact, print.level=0)
# ATT and its SE can be found in the est and se elements of m.exact
# Print the table
tab.exact = baltest.collect(mb.matched.exact, var.names = c('log_pop', 'age', 'educ_r', 'cashjob', 'urban', 'female', "home_afrik", "home_xhosa", "home_zulu") )
print(xtable(tab.exact[, 1:7], caption='Covariate Balance in Exactly Matched Data', label='tab:exact-bal'), caption.placement='top')
modelfits.t[5,] <- c(m.exact$est, m.exact$se, m.exact$orig.nobs, m.exact$wnobs)


##### START FIGURE 3 HERE ####
pdf(paste(fig, file="matchplots.pdf", sep=""), width = 9, height = 4)
par(mfrow = c(1,2), cex = .7)

plot(c(1,2,3,4,5), c(modelfits.w[1,1], modelfits.w[2,1],modelfits.w[3,1],modelfits.w[4,1],modelfits.w[5,1]),
     ylim = c(-.3,.3),
     xlim= c(0,6),
     xaxt = "n",
     xlab= "",
     ylab = "",
     pch = 16,
     main = "Water (Household)"
)

axis(tick = FALSE, 1, 1:5, labels = dvarlabels, cex.axis=1)

for (i in 1:5){
  lines(c(i,i), c(modelfits.w[i,1] - (1.96 * modelfits.w[i,2]),
                  modelfits.w[i,1] + (1.96 * modelfits.w[i,2])))
}

lines(c(0.0, 6.0), c(0,0), lty = 2)

plot(c(1,2,3,4,5), c(modelfits.t[1,1], modelfits.t[2,1],modelfits.t[3,1],modelfits.t[4,1],modelfits.t[5,1]),
     ylim = c(-.3,.3),
     xlim= c(0,6),
     xaxt = "n",
     xlab= "",
     ylab = "",
     pch = 16,
     main = "Toilet (Household/Compound)"
)

axis(tick = FALSE, 1, 1:5, labels = dvarlabels, cex.axis=1)

for (i in 1:5){
  lines(c(i,i), c(modelfits.t[i,1] - (1.96 * modelfits.t[i,2]),
                  modelfits.t[i,1] + (1.96 * modelfits.t[i,2])))
}

lines(c(0.0, 6.0), c(0,0), lty = 2)

dev.off()